Iskander et al, Science Translational Medicine, 2021

Figure S4 Characterization of EP phenotype and function.

load the library and object

library(SingCellaR)
library(ggplot2)
library(gridExtra)
library(ggpubr)

load(file = "/Users/gwang/DBA/DBA_revision/DBA_MKE.merged.louvain.rdata")

A) Expression of individual genes superimposed on FDG of all EP/MKP cells (n=6380) is shown, depicting a pattern of gene expression in normal mega-erythroid differentiation. Gray cells do not express the gene of interest. Color correlates with logarithmic-transformed expression as shown in color key.

plot_forceDirectedGraph_label_by_clusters(DBA_MKE,show_method = "merged_louvain")

plot_forceDirectedGraph_label_by_genes(DBA_MKE,gene_list = "CSF2RB")

plot_forceDirectedGraph_label_by_genes(DBA_MKE,gene_list = "CD34")

plot_forceDirectedGraph_label_by_genes(DBA_MKE,gene_list = "GATA1")

plot_forceDirectedGraph_label_by_genes(DBA_MKE,gene_list = "GATA2")

plot_forceDirectedGraph_label_by_genes(DBA_MKE,gene_list = "KLF1")

plot_forceDirectedGraph_label_by_genes(DBA_MKE,gene_list = "CD36")

plot_forceDirectedGraph_label_by_genes(DBA_MKE,gene_list = "ENG")

plot_forceDirectedGraph_label_by_genes(DBA_MKE,gene_list = "AHSP")

plot_forceDirectedGraph_label_by_genes(DBA_MKE,gene_list = "HBA1")

plot_forceDirectedGraph_label_by_genes(DBA_MKE,gene_list = "HBA2")

plot_forceDirectedGraph_label_by_genes(DBA_MKE,gene_list = "HBB")

plot_forceDirectedGraph_label_by_genes(DBA_MKE,gene_list = "EPOR")

plot_forceDirectedGraph_label_by_genes(DBA_MKE,gene_list = "PLEK")

plot_forceDirectedGraph_label_by_genes(DBA_MKE,gene_list = "PF4")

plot_forceDirectedGraph_label_by_genes(DBA_MKE,gene_list = "FLI1")

plot_forceDirectedGraph_label_by_genes(DBA_MKE,gene_list = "GP9")

plot_forceDirectedGraph_label_by_genes(DBA_MKE,gene_list = "VWF")

B) Violin plots show expression of E/MK-related transcription factors in control EP/MKP subclusters (n=2870; 1050 EEP, 393 Ecyc, 479 LEP, 948 MKP). Yellow dot indicates mean expression and the fraction of cells expressing each gene is depicted on the x axis.

plot_violin_for_genes_per_custom_group_of_cells<-function(object,custom_group_of_cells=list(),gene_list=c(),
                                                          take_log2=T,xlab.text.size=5,point.size=0.2,point.alpha=0.1,
                                                          grid.ncol = 1,grid.nrow=1){
  
  if(!is(object,"SingCellaR")){
    stop("Need to initialize the SingCellaR object")
  }
  if(length(gene_list)==0){
    stop("Required list of genes!")
  }
  
  ####################################
  umi.dat<-get_normalized_umi(object)
  my.p<-list()
  for(i in 1:length(gene_list)){
    my.gene<-as.character(gene_list)[i]
    my.sub.dat<-umi.dat[rownames(umi.dat)==my.gene,]
    my.new.dat<-data.frame()
    for(k in 1:length(custom_group_of_cells)){
      x.cell.names<-custom_group_of_cells[[k]]
      x.expr<-my.sub.dat[names(my.sub.dat) %in% as.character(x.cell.names)]
      ##################
      exp.val<-as.numeric(x.expr)
      n.exp<-length(exp.val[exp.val>0])
      n.total<-length(exp.val)
      my.label<-paste(names(custom_group_of_cells)[k],"\n",n.exp,"/",n.total,sep="")
      ##################
      y.legend<-""
      ##################
      if(take_log2==T){
        my.f<-data.frame(cluster=my.label,normUMI=log2(as.numeric(x.expr)+1))
        y.legend<-"log2(normalized UMI)"
      }else{
        my.f<-data.frame(cluster=my.label,normUMI=as.numeric(x.expr))
        y.legend<-"normalized UMI"
      }
      my.new.dat<-rbind(my.new.dat,my.f)
    }
    my.p[[i]]<-ggplot(
      data = my.new.dat,aes(x=cluster, y=normUMI,fill=cluster))+labs(y=y.legend)+
      geom_violin(trim=TRUE,scale="width")+
      geom_jitter(height = 0,size=point.size,alpha = point.alpha)+
      ggtitle(my.gene)+
      guides(fill=FALSE)+
      xlab("") +
      scale_fill_manual(values = c("#F8766D","#C77CFF","#00BFC4","#7CAE00")) +
      stat_summary(fun.y=mean, geom="point", size=1, color="yellow") +
      theme(axis.text=element_text(size=xlab.text.size))+
      theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank(),panel.background = element_blank(), axis.line = element_line(colour = "black"))
  } 
  grid.arrange(grobs = my.p,nrow=grid.nrow,ncol=grid.ncol)
}

umap.results <- DBA_MKE@umap.result
WT.cells <- umap.results$Cell[umap.results$genotype == "WT"]

sc.clusters <- DBA_MKE@sc.clusters
sc.clusters <- sc.clusters[sc.clusters$Cell %in% WT.cells,]
table(sc.clusters$merged_louvain_color,sc.clusters$merged_louvain)

##          
##            cl1  cl2  cl3  cl4
##   #00BFC4    0    0  479    0
##   #7CAE00    0  948    0    0
##   #C77CFF    0    0    0  393
##   #F8766D 1050    0    0    0

EEP <- sc.clusters$Cell[sc.clusters$merged_louvain == "cl1"]
LEP <- sc.clusters$Cell[sc.clusters$merged_louvain == "cl3"]
Ecyc <- sc.clusters$Cell[sc.clusters$merged_louvain == "cl4"]
MKP <- sc.clusters$Cell[sc.clusters$merged_louvain == "cl2"]



plot_violin_for_genes_per_custom_group_of_cells(DBA_MKE,
                                                custom_group_of_cells = list(EEP,Ecyc,LEP,MKP),
                                                gene_list = "GATA1",
                                                take_log2 = T
)

plot_violin_for_genes_per_custom_group_of_cells(DBA_MKE,
                                                custom_group_of_cells = list(EEP,Ecyc,LEP,MKP),
                                                gene_list = "GATA2",
                                                take_log2 = T
)

plot_violin_for_genes_per_custom_group_of_cells(DBA_MKE,
                                                custom_group_of_cells = list(EEP,Ecyc,LEP,MKP),
                                                gene_list = "KLF1",
                                                take_log2 = T
)

plot_violin_for_genes_per_custom_group_of_cells(DBA_MKE,
                                                custom_group_of_cells = list(EEP,Ecyc,LEP,MKP),
                                                gene_list = "FLI1",
                                                take_log2 = T
)